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Abstract 

We briefly review advances in understanding the initial stages of a heavy ion collision. In particular the focus is on 
moving from parametrizing the initial state to calculating its properties from QCD, consistently with the description 
of hard probes and dilute-dense scattering experiments. Modeling the event-by-event fluctuating nuclear geometry in 
initial state calculations has significantly improved in recent years. We also discuss prospects of directly seeing effects 
of particle correlations created in the initial state in the experimental observables. 
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1. Introduction 

For a practitioner of hydrodynamic modeling of 
heavy ion collisions, the question of initial conditions 
often boils down to an acronym soup of models. One 
tries out a suitable subset of these (“Glauber”, KLN, 
[mc]KLN, mcrcBK, EPOS, EKRT, IPglasma etc.) as 
an input to a fluid dynamics calculation and compares 
to experimental data. The purpose here is very differ¬ 
ent, namely to discuss the physics ingredients in these 
calculations and concentrate more on their common as¬ 
pects them than on the differences. This is done with 
the caveat that we are concentrating exclusively on a 
weak coupling partonic description of the degrees of 
freedom involved. We will start with a very brief in¬ 
troduction to the physics picture in the weak coupling 
approach. We then discuss dilute-dense control exper¬ 
iments that can be used to more directly probe the de¬ 
grees of freedom in the initial state of a heavy ion col¬ 
lision. The most recent dynamical initial state models 
will then be discussed in comparison to Monte Carlo 
Glauber parametrizations. Einally we will describe re¬ 
cent calculations of long range correlations that, if one 
is able to discriminate between effects of collective flow 
and those originating from the initial state, could pro¬ 
vide a direct access to the latter. 



Figure 1: Gluon cascade leading to particle production in the central 
rapidity region. 

2. Initial state at weak coupling 

The weak coupling picture of particle production at 
central rapidities in very high energy hadronic collisions 
is depicted in Eig. 1 . The gluons that are mainly respon¬ 
sible for particle production result from a cascade of 
multiple splittings from the initial valence quarks. The 
emission probability for one splitting is q-s dx/.r, where 
X is the fraction of the longitudinal momentum of the in¬ 
coming hadron and as is the QCD coupling. This form 
is constant as a function of rapidity, and thus leads nat¬ 
urally to a rapidity plateau, i.e. a multiplicity dV/ dy 
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approximately independent of y for scales Ay < I/cKs- 
While at RHIC the total collision energy is still so low 
that the y-dependence is dominated by large-;*: physics, 
at the LHC this plateau has become prominently visible 
in the experimental data. At high enough energy there 
is phase space for many gluon emissions in the cascade, 
each additional one suppressed by a factor but en¬ 
hanced by a phase space volume Ay ~ In ^fs. With a 
factor 1/n! from the rapidity ordering of n gluons one 
could in fact very roughly expect 
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gluons. Again, the experimental energy dependence of 
multiplicities in pp (~ and AA (~ colli¬ 

sions fits in very nicely with this very crude weak cou¬ 
pling, small angle scattering estimate. This could be 
contrasted with a strong coupling picture, which would 
genetically lead to complete baryon stopping [1] and a 
stronger growth of the multiplicity with ^fs. 

Eventually the gluons in this cascade will overlap and 
gluon mergings will start to compete with splittings. 
This increase in the gluon density corresponds to the 
Yang-Mills theory becoming completely nonperturba- 
tive. This happens when the two terms of the covari¬ 
ant derivative -iD/j = -idf^ + gA^ - + gA^ are of 

the same order. The (transverse) momentum scale at 
which this happens is referred to as the saturation scale 
Pt ~ gAx ~ Qs- When the energy is high enough, 
2s ^ Aqcd and the coupling is still weak. What re¬ 
sults is a picture that has a weak coupling g but is still 
nonperturbative due to the large gluon field strengths 
~ !/.?■ Since the occupation numbers of gluonic 
states f{k) ~ ~ 1 ja^ are large, the gluon field is, 

to leading order in the coupling, a classical field. 

Our understanding of how exactly such a system of 
strong (but very anisotropic) gluon fields thermalizes (or 
slightly more modestly isotropizes) has seen significant 
progress in the recent years (for a review see [2]). Pur¬ 
suing this question in any depth is out of scope here, 
but as a rough outline the process must behave in the 
following way. The initial classical fields, with with oc¬ 
cupation numbers f{k) ~ ^ j ^re very anisotropic 

due to the longitudinal expansion of the system. As glu¬ 
ons start to scatter off the transverse plane and occupy 
a larger volume in phase space, the occupation number 
decreases. When f{k) <si ^ one can switch to a kinetic 
theory description and follow the system towards local 
equilibrium. Close to isotropy kinetic theory matches 
smoothly into a order viscous hydrodynamical de¬ 
scription. A recent weak coupling calculation using an 


Figure 2: Different fragmentation functions in NLO pQCD compared 
to CMS data for charged particle spectra in proton-proton collisions, 
from [9]. 


effective kinetic theory for QCD [3] arrives, when ex¬ 
trapolating to realistic values of as, to times as short as 
1 fm for a matching to viscous hydrodynamics. 


3. Dilute-dense control measurements 

The quintessential experiment for measuring the par- 
tonic content of a hadron or nucleus is deep inelastic 
scattering (DIS). A good description of the initial state 
of a heavy ion collision should be consistent with the 
precise measurements of quarks and gluons in a proton 
at HERA and with the available electron-nucleus cross 
section data. It should also be able to give quantita¬ 
tive predictions for future measurements at an EIC. This 
is naturally true for calculations whose starting point 
are nuclear pdfs (e.g. EPS09 [4]). It also applies to 
more recent CGC calculations using e.g. the IPSat or 
bCGC parametrizations [5, 6] and to calculations us¬ 
ing the mnning coupling Balitsky-Kovchegov equation 
(rcBK), e.g. [7, 8]. In the classical field CGC picture the 
initial color fields in a heavy collision are calculated in 
terms of Wilson lines, whose correlator determines the 
total DIS cross section. 

Another check of QCD descriptions is provided by 
hadronic dilute-dense, i.e. proton-nucleus, collisions. 
The most straightforward observables here are ratios of 
particle spectra in proton-nucleus collisions to those in 
proton-proton ones, normalized with the number of nu¬ 
cleons in the nucleus, known as RpA- At high pr these 
ratios would be expected to approach unity for many 
particle species. This is particularly important when 
moving from minimum bias proton-nucleus collisions 
to separate centrality classes, where hard piticle pro¬ 
duction serves as an important consistency check of the 
Glauber modeling [10] needed for the centrality deter¬ 
mination. If the normalization of Rp^ is under control. 
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its small-/? 7 - behavior is sensitive to the enhancement 
of saturation effects in a nucleus, especially at forward 
rapidity where one probes smaller values of x in the nu¬ 
cleus. 

Typical CGC calculations of hadronic RpA at semi- 
hard pt use two different formalisms. The forward 
rapidity “hybrid” picture starts from a collinear quark 
or gluon from the probe, which is described in terms 
of conventional DGLAP evolved parton distribution 
functions. This parton propagates through the dense 
color field of the target, picking up an eikonal Wilson 
line. The other possible treatment is a symmetric kr- 
factorized description, where the produced gluon spec- 
turn is obtained as a convolution of two unintegrated 
gluon distributions (ugd’s) describing the projectile and 
the target. The ugd is, again, obtained as a Fourier- 
transform of the Wilson line correlator probed in DIS. 
Both formalisms require, unexceptionally for leading 
order QCD calculations, a “/T-factor” of order ~ 2 to fit 
the spectra, but this normalization uncertainty cancels 
in the ratio RpA- Unlike in the early days of RHIC, up 
to date calculations of RpA [11, 12, 13, 8] use ugd’s that 
have been fit to the precise HERA data. In fact, because 
this data is very constraining, the main quantitative dif¬ 
ference between these calculations is the treatment of 
the nuclear geometry needed to go from a proton to a 
nuclear target. With RpA predictions differing, e.g. by 
a factor ~ 2 from Ref. [11] to [8], this can be a signif¬ 
icant effect. A similar example of the importance of a 
proper treatment of the nuclear geometry is provided by 
forward J/'¥ production [14]. 

Even more discriminative power is in principle pro¬ 
vided by the absolute spectra themselves. In a weak 
coupling QCD calculation the normalization has, how¬ 
ever, much more uncertainty than the ratio RpA . An ad¬ 
ditional complication comes from the recent observa¬ 
tion [9] that current fragmentation functions do not pro¬ 
vide a very good description of data (see Eig. 2) even 
in the cases where they should work well, i.e. proton- 
proton collisions at higher pr- This will require fur¬ 
ther development, focusing in particular on the poorly 
constrained gluon fragmentation funtions, see e.g. [15]. 
Electroweak probes (photons and electroweak bosons) 
are theoretically the cleanest observable probing weak 
coupling partonic physics in the nucleus. Here the 
greater question is whether the experimental accuracy 
will be good enough to distinguish nuclear effects. 

4. MC Glauber and dynamical models 

Initial matter formation in heavy ion collisions is of¬ 
ten parametrized in terms of a “Monte Carlo Glauber” 


(MCG) model [16]. One starts from nucleons dis¬ 
tributed inside the nucleus according to a standard 
Woods-Saxon nuclear density profile. Nucleons from 
the colliding projectiles are then deemed to “collide” 
or “participate” if they lie within a short enough dis¬ 
tance (determined by the total inelastic nucleon-nucleon 
cross section) from each other. Experimental event-by 
event distributions of charged particle multiplicities or 
forward calorimeter energies can be roughly modeled 
by assuming that the multiplicity or energy is propor¬ 
tional to the number of “participant nucleons” or “bi¬ 
nary nucleon-nucleon collisions”. This has led to the 
MCG model being used both experimentally in event 
centrality classification, and theoretically as a spatial 
distribution of the initial conditions of hydrodynamical 
simulations. 

In hydrodynamical calculations microscopical mod¬ 
els of particle production are sometimes contrasted with 
MCG based on comparisons to experimental data. We 
would like to argue that this is a fundamentally mis¬ 
leading and unfair comparison. In addition to not be¬ 
ing completely unique (one can e.g. choose the mul¬ 
tiplicity or the energy as proportional to either Acoii or 
A^part), the MCG model is a purely empirical observation 
without any dynamical mechanism for particle produc¬ 
tion. It does not, for example, make any predictions 
about rapidity or dependence. The MCG model is 
such a successful description of nuclear geometry that 
a similar picture is in fact built into most modern initial 
state calculations which combine it with a microscopi¬ 
cal description of particle production (strings, classical 
fields, parton scattering ...). Thus it is not surprising 
that, to a first approximation, the centrality dependence 
of bulk observables from these models is close to MCG. 
On a conceptual level, if one is interested in understand¬ 
ing microscopical QCD dynamics, the meaningful com¬ 
parison is between different dynamical models of parti¬ 
cle production, not between them and the MCG. One 
can ask whether the theory is consistent within itself, 
or how well the (necessarily small) deviations from a 
dNch/ dr] ~ Apart scaling that it predicts describe experi¬ 
mental observations. 

Let us recall some examples of how a very simi¬ 
lar treatment of the nuclear geometry is at the heart 
of many recent initial state calculations. The MCKLN 
model [17] explicitly starts from an MCG calculation 
that determines, event-by-event, a set of participant nu¬ 
cleons in each nucleus. One then takes the (now probe- 
dependent, i.e. nonuniversal) saturation scales in the nu¬ 
clei as proportional to this Apart and uses them as inputs 
in a kj-factorized calculation of the initial energy den¬ 
sity. The MCrcBK model [18] is otherwise similar, but 
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the ad hoc functional form of the ugd is replaced by a 
solution of the BK equation. The IPglasma initial state 
model [19] starts from the IPsat [20] parametrization for 
the dipole scattering amplitude in a nucleus 

= I _ (^2) 

involving fluctuating nucleon positions bn from a 
Woods-Saxon distribution. One then determines a local 
saturation scale 2s(l>r) from this amplitude and uses it 
in an MV-model correlator of color charges 

<p“(x7-)p*(yr)) ~ 2s(xr)5“''d®(x7- - yr), (3) 

from which the classical Yang-Mills (CYM) field is cal¬ 
culated. The EKRT model in its present version [21] 
uses a fluctuating nuclear thickness profile 

A 

TA(bT)^Yj^pibT-bn) (4) 

1=1 

to determine the initial energy density as 

e(bT) = r[TA{bT)TB{bT)], (5) 


use a ^ 7 --factorized formula to compute gluon produc¬ 
tion the procedure is slightly different. Because of par- 
ton saturation the total gluon multiplicity is only loga¬ 
rithmically infrared divergent 

—-r =- 2 f ‘Py(^T)‘Py(PT - ^r), (7) 

d^pj. dy as pP Jkr 

where i^y(k 7 -) is the unintegrated gluon distribution. 
Thus, in stead of determining the cutoff together with 
the normalization, one usually chooses an infrared cut¬ 
off scale (often by restricting the kj-integral to kr < pn 
for a different scheme see [12]) and then adjusts the nor¬ 
malization constant C to data. This leaves the gluon 
mean pj (energy per particle) quite unconstrained and 
potentially very far from a value extrapolated backward 
from the experimental data [18]. 

In the CYM calculations, such as the IPglasma 
model [19], the gluon multiplicity is finite due to nonlin¬ 
ear effects in the final state [25] and no infrared cutoffs 
are needed. The Wilson lines U{xt) that fully deter¬ 
mine the initial conditions of the CYM calculation [26] 
are constrained by DIS data because the correlator 


where the QCD dynamics in the calculation is 
parametrized by the function !F. Also string-based mod¬ 
els such as EPOS [22] and AMPT [23] share the same 
description of the nuclear geometry. Since we are still 
far from being able to calculate nuclear structure from 
first principles QCD, at some level the MCG geometry 
has to be put in by hand, no matter the sophistication of 
the QCD calculation 

Besides the spatial, distribution the calculation of the 
initial state should naturally also provide the value of 
the energy (or entropy) density. Ideally this should be 
done without any adjustable free parameters, but as we 
will now discuss this is not fully the case for any of the 
currently available calculations. In a purely perturbative 
calculation the initial gluon multiplicity is not finite, so 
the question of normalization is intimately tied in with 
the issue of infrared cutoffs. This is a particularly salient 
feature of the EKRT model [24], where the initial parton 
spectrum calculated in collinear perturbation theory has 
a strong power law dependence on the infrared cutoff. 
The model introduces two essential free parameters, the 
cutoff (saturation scale) po and a normalization coeffi¬ 
cient Ksat, which are simultaneously determined from 
the hnal state saturation criterion 


dEripr > Po) 



n 


( 6 ) 


and by fitting the multiplicity in the most central LHC 
lead-lead collisions. In KLN-type initial conditions that 


-^TrUHxTWiyr) ( 8 ) 

is proportional to the total DIS cross section. Thus on 
a conceptual level there are no free adjustable parame¬ 
ters in the CYM scenario. The practical implementation 
in IPglasma, however, uses a procedure that leaves room 
for small uncertainties. One first uses a fit to HERA data 
to get the saturation scale Q^ibr)- This is then input 
into an MV model calculation so that one reproduces 
the same Qs, but not the full Wilson line correlator of 
the DIS ht. Additionally, in the MV model in a finite 
nucleus one must regulate long distance Coulomb tails 
of the classical gauge helds with a confinement scale 
parameter. The net result of this procedure is that the 
Wilson line in the CYM initial condition does not ex¬ 
actly match the original DIS cross section, although it is 
very well constrained by it. 


5. Correlations: direct signals of the initial state? 

A simple causality argument shows that correlations 
in particle production that have a long range in rapidity 
must have originated in the earliest stage of the collision 
process. This is analogous to the way correlations in 
the cosmic microwave background are sensitive to early 
times, when now separate parts of the universe were still 
in causal contact. Azimuthal flow coefficients v„ are an 
important example of long range rapidity correlations. 
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Figure 3: Color field domains in the target 

They are either explicitly constructed from a (usually ra¬ 
pidity separated) multiparticle correlation function (the 
“cumulant method”) or by correlating particles at one 
rapidity with a reaction plane determined from particles 
in another part of the detector. The same correlation can 
be presented in terms of a yield per trigger or as a flow 
coefficient v„; see e.g. [27]. The most important origin 
of these correlations, especially in large collision sys¬ 
tems, is the geometry of the initial state; the distribution 
of matter in the transverse plane is by definition felt at 
all rapidities. Geometrical correlations are, however, in 
position space, and strong enough collective final state 
interactions are needed to transform them into observ¬ 
able momentum space distributions. 

It has more recently [28, 29] been pointed out that, in 
particular in small collision systems, the domain struc¬ 
ture in the target color field can also generate azimuthal 
correlations directly in momentum space without final 
state collective phenomena. These correlations follow 
in a very intuitive way from the CGC picture. The usual 
picture of particle production in this case is that of an 
individual quark or gluon (whose number is given by a 
conventional pdf) passing through the strong color field 
of the target and being deflected to some transverse mo¬ 
mentum. Since the target color field consists of domains 
of size ~ 1 /2s (s®® Fig- 3), incoming particles in the 
same color state and hitting the same domain experi¬ 
ence a similar deflecting color field and become corre¬ 
lated. The mechanism generates correlations that are 
suppressed by the number of independent domains and 
colors as ~ ^liN^Q^S j_), where 5^ is the size of inter¬ 
action area. Thus, in contrast to collective flow effects, 
these correlations are enhanced in small collision sys¬ 
tems. 

These correlations have recently been analyzed by 
different authors in calculations that share the same 
physical picture, but differ in the approximations used 


to calculate the correlations of the target color fields. 
The “Glasma graph” calculations [28, 30] linearize in 
the color charge density of the target and assume Gaus¬ 
sian correlations between the domains. In the “color 
field domain model” [31] one also linearizes in the tar¬ 
get field, but adds an additional non-Gaussian correla¬ 
tion. On the other hand one can perform a full nonlin¬ 
ear calculation in the dilute-dense limit [32, 33] with¬ 
out adding non-Gaussian correlations. By performing 
a full CYM simulation [34] one can also include pre¬ 
equilibrium collective effects in the calculation of the 
azimuthal dependence. For a more detailed reecnt dis¬ 
cussion of the differences and similarities see [35]. So 
far all of these calculations assume parametrically small 
rapidity separations Ay < I/a^, but there have been pro¬ 
posals to extend the calculation also to parametrically 
large rapidity intervals [36]. 

6. Conclusions 

In conclusion, recent years have seen a huge progress 
towards a more well defined weak coupling picture of 
the initial stages of a heavy ion collision. Studies of 
thermalization are moving from qualitative to quanti¬ 
tative, and the importance of interfacing with kinetic 
theory for the later stage of equilibration is becoming 
understood. One wants the initial stage calculation to 
be consistent with perturbative probes and control mea¬ 
surements in dilute-dense collisions. We emphasized 
that the MC Glauber-like geometry is built into the more 
modern dynamical models of the initial state. Finally 
we discussed the azimuthal structure of long range ra¬ 
pidity correlations, i.e. the v„-coefficients, especially in 
small systems, where the interplay between initial and 
final state collective effects remains to be quantitatively 
sorted out. 

This work has been supported by the Academy of 
Finland, projects 267321 and 273464. 
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